rm(list=ls())
gc()

## ESTIMATE TSMART PID PLACEBO PRE-PERIOD MODELS

## LOAD PACKAGES
library(data.table)
library(lfe)
library(stringr)



data = fread('panel-analysis-2016-2020.csv.gz',select = c('DemSpExpDiff_nohh', 'RepSpExpDiff_nohh','DemSpExp_nohh_year1','RepSpExp_nohh_year1','Age_year1','Race','Gender',
                                                                           'ZipCode','countyfips','Party_year1','hh.n.adj_year1','hh.n.adj_year2',
                                                                           'hh.d.adj_year1','hh.d.adj_year2','hh.r.adj_year1','hh.r.adj_year2',
                                                                           'Married_year1','State', 'WhiteBlockGroupDiff','MarriageDiff','AgeBlockGroupDiff','RegsBlockGroupDiff',
                                                                           'HHIncomeBlockGroupDiff' , 'CollegeBlockGroupDiff' , 'HomeownerBlockGroupDiff' , 'YearBuiltBlockGroupDiff' , 
                                                                           'DriveWorkBlockGroupDiff' , 'EmplBlockGroupDiff' , 'HouseValueBlockGroupDiff', 'Party_2015','Party_2014','Party_2013','Party_2012'))

data[!Party_year1%in%c('Democrat','Republican'), Party_year1:='Other']
data[countyfips=='',countyfips:=NA]
data[ZipCode=='',ZipCode:=NA]
data[Gender=='',Gender:=NA]

# Make household composition variables
data[,hh.n.diff:=hh.n.adj_year2-hh.n.adj_year1]
data[,hh.d.diff:=hh.d.adj_year2-hh.d.adj_year1]
data[,hh.r.diff:=hh.r.adj_year2-hh.r.adj_year1]


# Empty pre-period party values (i.e., voters for whom we have no information on that year's past partisanship) should be coded as NAs so they are dropped from the estimation for that year
data[Party_2012=='',Party_2012:=NA]
data[Party_2013=='',Party_2013:=NA]
data[Party_2014=='',Party_2014:=NA]
data[Party_2015=='',Party_2015:=NA]

# Make pre-period party variables
data[!is.na(Party_2012)&!Party_2012%in%c('Democrat','Republican'),Party_2012:='Non-Partisan']
data[!is.na(Party_2013)&!Party_2013%in%c('Democrat','Republican'),Party_2013:='Non-Partisan']
data[!is.na(Party_2014)&!Party_2014%in%c('Democrat','Republican'),Party_2014:='Non-Partisan']
data[!is.na(Party_2015)&!Party_2015%in%c('Democrat','Republican'),Party_2014:='Non-Partisan']


data[,Dem_2012 := as.numeric(Party_2012=='Democrat')]
data[,Dem_2013 := as.numeric(Party_2013=='Democrat')]
data[,Dem_2014 := as.numeric(Party_2014=='Democrat')]
data[,Dem_2015 := as.numeric(Party_2015=='Democrat')]


data[,Rep_2012 := as.numeric(Party_2012=='Republican')]
data[,Rep_2013 := as.numeric(Party_2013=='Republican')]
data[,Rep_2014 := as.numeric(Party_2014=='Republican')]
data[,Rep_2015 := as.numeric(Party_2015=='Republican')]

  # split data
  dems = data[Party_year1=='Democrat']
  reps = data[Party_year1=='Republican']
  oths = data[Party_year1=='Other']
  
  rm(data)
  gc()
  
  # construct grouping variables based on age decile, gender, race, year1 treatment decile, and city
      # we can drop files where any of these are missing

          
dems = dems[!is.na(Age_year1) & !is.na(Gender) & !is.na(Race) & !is.na(Married_year1) &!is.na(ZipCode)]
reps = reps[!is.na(Age_year1) & !is.na(Gender) & !is.na(Race) & !is.na(Married_year1) & !is.na(ZipCode)]
oths = oths[!is.na(Age_year1) & !is.na(Gender) & !is.na(Race) & !is.na(Married_year1) & !is.na(ZipCode)]

dems[,AgeDecile:=cut(Age_year1,breaks=unique(as.numeric(quantile(Age_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]
reps[,AgeDecile:=cut(Age_year1,breaks=unique(as.numeric(quantile(Age_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]
oths[,AgeDecile:=cut(Age_year1,breaks=unique(as.numeric(quantile(Age_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]


# spatial seg treatment
# democrat exposure
dems[,DemSpDecile:=cut(DemSpExp_nohh_year1,breaks=unique(as.numeric(quantile(DemSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]
reps[,DemSpDecile:=cut(DemSpExp_nohh_year1,breaks=unique(as.numeric(quantile(DemSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]
oths[,DemSpDecile:=cut(DemSpExp_nohh_year1,breaks=unique(as.numeric(quantile(DemSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]

# republican exposure
dems[,RepSpDecile:=cut(RepSpExp_nohh_year1,breaks=unique(as.numeric(quantile(RepSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]
reps[,RepSpDecile:=cut(RepSpExp_nohh_year1,breaks=unique(as.numeric(quantile(RepSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]
oths[,RepSpDecile:=cut(RepSpExp_nohh_year1,breaks=unique(as.numeric(quantile(RepSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))),include.lowest=T)]



#####
# make groups
# spatial seg treatment
dems[,GroupSpDemExp:=paste(ZipCode,AgeDecile,Race, Gender, DemSpDecile,State,hh.n.adj_year1,hh.d.adj_year1,Married_year1, sep ='_')]
reps[,GroupSpDemExp:=paste(ZipCode,AgeDecile,Race, Gender, DemSpDecile,State,hh.n.adj_year1,hh.d.adj_year1,Married_year1, sep ='_')]
oths[,GroupSpDemExp:=paste(ZipCode,AgeDecile,Race, Gender, DemSpDecile,State,hh.n.adj_year1,hh.d.adj_year1,Married_year1, sep ='_')]

dems[,GroupSpRepExp:=paste(ZipCode,AgeDecile,Race, Gender, RepSpDecile,State,hh.n.adj_year1,hh.r.adj_year1,Married_year1, sep ='_')]
reps[,GroupSpRepExp:=paste(ZipCode,AgeDecile,Race, Gender, RepSpDecile,State,hh.n.adj_year1,hh.r.adj_year1,Married_year1, sep ='_')]
oths[,GroupSpRepExp:=paste(ZipCode,AgeDecile,Race, Gender, RepSpDecile,State,hh.n.adj_year1,hh.r.adj_year1,Married_year1, sep ='_')]

###
  # estimate models
  
# 2012
ModelDemSpExpDems_2012 = summary(felm(Dem_2012 ~ DemSpExpDiff_nohh + hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpReps_2012 = summary(felm(Dem_2012 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpOths_2012 = summary(felm(Dem_2012 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)

ModelRepSpExpDems_2012 = summary(felm(Rep_2012 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpReps_2012 = summary(felm(Rep_2012 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpOths_2012 = summary(felm(Rep_2012 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)



#2013
ModelDemSpExpDems_2013 = summary(felm(Dem_2013 ~ DemSpExpDiff_nohh + hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpReps_2013 = summary(felm(Dem_2013 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpOths_2013 = summary(felm(Dem_2013 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)

ModelRepSpExpDems_2013 = summary(felm(Rep_2013 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpReps_2013 = summary(felm(Rep_2013 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpOths_2013 = summary(felm(Rep_2013 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)


#2014
ModelDemSpExpDems_2014 = summary(felm(Dem_2014 ~ DemSpExpDiff_nohh + hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpReps_2014 = summary(felm(Dem_2014 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpOths_2014 = summary(felm(Dem_2014 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)

ModelRepSpExpDems_2014 = summary(felm(Rep_2014 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpReps_2014 = summary(felm(Rep_2014 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpOths_2014 = summary(felm(Rep_2014 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)


#2015
ModelDemSpExpDems_2015 = summary(felm(Dem_2015 ~ DemSpExpDiff_nohh + hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpReps_2015 = summary(felm(Dem_2015 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpOths_2015 = summary(felm(Dem_2015 ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)

ModelRepSpExpDems_2015 = summary(felm(Rep_2015 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpReps_2015 = summary(felm(Rep_2015 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpOths_2015 = summary(felm(Rep_2015 ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+MarriageDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                        HHIncomeBlockGroupDiff + CollegeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                        DriveWorkBlockGroupDiff + EmplBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)



#

save(ModelDemSpExpDems_2015, ModelDemSpExpReps_2015, ModelDemSpExpOths_2015, 
     ModelRepSpExpDems_2015, ModelRepSpExpReps_2015, ModelRepSpExpOths_2015,
     
     ModelDemSpExpDems_2014, ModelDemSpExpReps_2014, ModelDemSpExpOths_2014, 
     ModelRepSpExpDems_2014, ModelRepSpExpReps_2014, ModelRepSpExpOths_2014,
     
     ModelDemSpExpDems_2013, ModelDemSpExpReps_2013, ModelDemSpExpOths_2013, 
     ModelRepSpExpDems_2013, ModelRepSpExpReps_2013, ModelRepSpExpOths_2013,
     
     ModelDemSpExpDems_2012, ModelDemSpExpReps_2012, ModelDemSpExpOths_2012, 
     ModelRepSpExpDems_2012, ModelRepSpExpReps_2012, ModelRepSpExpOths_2012,
     
     file = 'results/placebo-results-2016-2020.Rdata')

